Activated and renewed HBC clustering and comparison with scATAC-seq
Koen Van den Berge
7/22/2020
here::set_here()## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scriptssuppressPackageStartupMessages({
library(slingshot)
library(SingleCellExperiment)
library(cowplot)
library(rgl)
library(clusterExperiment)
library(RColorBrewer)
library(ggplot2)
library(pheatmap)
library(wesanderson)
library(gridExtra)
library(irlba)
library(tidyverse)
})Import data
sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
# counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
load("../data/ALL_TF.Rda")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")
timePoint <- readRDS("../data/finalTrajectory/timePoint_noResp_noMV.rds")Select activated and renewed HBC clusters
I will select
- activated HBCs as cells sampled in the 24h timepoint that are in the starting cluster
- renewed (‘resting’) HBCs as cells sampled in the 14D timepoint that are in the HBC cluster
cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "#999999")
view3d( theta = 10, phi = 10)
# Datta labels
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clDatta], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")
# timepoints
plot3d(reducedDim(sds), aspect = 'iso', col=cols[timePoint], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")
# cluster labels used for TI
clSds <- apply(slingClusterLabels(sds), 1, which.max)
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clSds], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")
actHBCId <- which(clSds == 8 & timePoint == "24HPI")
restHBCId <- which(clDatta == "HBC" & timePoint == "14DPI")
grpHBC <- rep(NA, length(c(actHBCId, restHBCId)))
grpHBC[1:length(actHBCId)] <- "activated"
grpHBC[is.na(grpHBC)] <- "resting"
grpHBC <- factor(grpHBC)
countHBC <- counts[, c(actHBCId, restHBCId)]Dimensionality reduction
library(scran)
library(uwot)
library(scater)##
## Attaching package: 'scater'## The following object is masked from 'package:clusterExperiment':
##
## plotHeatmaplogFracs <- log(sweep(countHBC+1e-5, 2, colSums(countHBC), "/"))
gv <- modelGeneVar(countHBC)
hvg <- getTopHVGs(gv, n=1000)
pc10 <- irlba::irlba(logFracs[hvg,], nv=10)
um10 <- uwot::umap(pc10$v %*% diag(pc10$d))
plot(um10, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")
Batch correction with Harmony
library(harmony)## Loading required package: Rcppharmony_embeddings <- HarmonyMatrix(logFracs, meta_data=grpHBC, do_pca=TRUE)## Harmony 1/10## Harmony 2/10## Harmony converged after 2 iterationsumHarmony <- uwot::umap(harmony_embeddings)
plot(umHarmony, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")
Integration: Processing and dimensionality reduction of each dataset separately
library(Signac)
library(Seurat)##
## Attaching package: 'Seurat'## The following object is masked from 'package:SummarizedExperiment':
##
## Assayslibrary(patchwork)##
## Attaching package: 'patchwork'## The following object is masked from 'package:cowplot':
##
## align_plotslibrary(ArchR)## Loading required package: data.table##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transpose## The following object is masked from 'package:SummarizedExperiment':
##
## shift## The following object is masked from 'package:GenomicRanges':
##
## shift## The following object is masked from 'package:IRanges':
##
## shift## The following objects are masked from 'package:S4Vectors':
##
## first, second## Loading required package: rhdf5## Loading required package: magrittr##
## Attaching package: 'magrittr'## The following object is masked from 'package:purrr':
##
## set_names## The following object is masked from 'package:tidyr':
##
## extractset.seed(1234)
addArchRGenome("mm10")## Setting default genome to Mm10.# scATAC-seq importing and preprocessing
# from 20200610_scATAC_injury_archr_samples_pairwiseHBC_v3.Rmd
archProj <- loadArchRProject("../data/scATAC/archR_samples_pairwise_v2")## Successfully loaded ArchRProject!##
## / |
## / \
## . / |.
## \\\ / |.
## \\\ / `|.
## \\\ / |.
## \ / |\
## \\#####\ / ||
## ==###########> / ||
## \\##==......\ / ||
## ______ = =|__ /__ || \\\
## ,--' ,----`-,__ ___/' --,-`-===================##========>
## \ ' ##_______ _____ ,--,__,=##,__ ///
## , __== ___,-,__,--'#' ===' `-' | ##,-/
## -,____,---' \\####\\________________,--\\_##,/
## ___ .______ ______ __ __ .______
## / \ | _ \ / || | | | | _ \
## / ^ \ | |_) | | ,----'| |__| | | |_) |
## / /_\ \ | / | | | __ | | /
## / _____ \ | |\ \\___ | `----.| | | | | |\ \\___.
## /__/ \__\ | _| `._____| \______||__| |__| | _| `._____|
## clHBCMerged <- archProj$clHBCMerged
plotEmbedding(ArchRProj = archProj, colorBy = "cellColData",
name = "clHBCMerged", embedding = "UMAP_cor")## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-175b8618f04a2-Date-2020-12-03_Time-17-43-46.log
## If there is an issue, please report to github with logFile!## Getting UMAP Embedding## ColorBy = cellColData## Plotting Embedding## 1
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-175b8618f04a2-Date-2020-12-03_Time-17-43-46.log
archrUmap <- ArchR::getEmbedding(archProj, embedding = "UMAP_cor")
colnames(archrUmap) <- c("UMAP1", "UMAP2")
archrUmap$clHbcMerged <- getCellColData(archProj, "clHBCMerged")$clHBCMerged
archrUmap$treatment <- getCellColData(archProj, "treat")$treat
peaks <- archProj@peakSet
names(peaks) <- paste0("peak",1:length(peaks))
atacPeakMat <- readRDS("../data/scATAC/peakMatArchrHbc_v2.rds")
rownames(atacPeakMat) <- names(peaks)
geneActivity <- readRDS("../data/scATAC/geneMatArchrHbc.rds")
atac <- CreateSeuratObject(counts = atacPeakMat,
assay = "ATAC",
project = "10x_ATAC")
atac[["ACTIVITY"]] <- CreateAssayObject(counts = geneActivity)
atac$tech <- "atac"
atac$treatment <- archProj@cellColData$treat
atac$clHBCMerged <- clHBCMerged
# process gene activity
DefaultAssay(atac) <- "ACTIVITY"
atac <- FindVariableFeatures(atac)
atac <- NormalizeData(atac)
atac <- ScaleData(atac)## Centering and scaling data matrix# process peak matrix
DefaultAssay(atac) <- "ATAC"
# VariableFeatures(atac) <- names(which(Matrix::rowSums(atac) > 50))
# atac <- RunLSI(atac, n = 30, scale.max = NULL)
# atac <- RunUMAP(atac, reduction = "lsi", dims = 1:30)
atac <- RunTFIDF(atac)## Performing TF-IDF normalizationatac <- FindTopFeatures(atac, min.cutoff = 'q0')
atac <- RunLSI(atac, n = 30)## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac
## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac
## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac## Performing TF-IDF normalization
## Running SVD on TF-IDF matrix
## Scaling cell embeddingsatac <- RunUMAP(object = atac, reduction = 'lsi', dims = 2:30)## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 17:45:07 UMAP embedding parameters a = 0.9922 b = 1.112
## 17:45:07 Read 4076 rows and found 29 numeric columns
## 17:45:07 Using Annoy for neighbor search, n_neighbors = 30
## 17:45:07 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 17:45:08 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpBeaM9s/file175b8747959b7
## 17:45:08 Searching Annoy index using 1 thread, search_k = 3000
## 17:45:09 Annoy recall = 100%
## 17:45:10 Commencing smooth kNN distance calibration using 1 thread
## 17:45:12 Initializing from normalized Laplacian + noise
## 17:45:12 Commencing optimization for 500 epochs, with 152690 positive edges
## 17:45:19 Optimization finished# scRNA-seq import
rna <- CreateSeuratObject(counts=countHBC)
rna$tech <- '10x'
rna$celltype <- grpHBC
rna <- FindVariableFeatures(rna, selection.method = "vst", nfeatures = 2000)
rna <- ScaleData(rna)## Centering and scaling data matrixrna <- RunPCA(rna, npcs=20)## PC_ 1
## Positive: Slc6a6, Rps27, Junb, Aqp3, Btg2, Fos, Fxyd3, Apoe, Krt5, Jun
## Gclm, Socs3, Ier2, Gpm6a, Ubc, Ctsl, Cyr61, Ifitm3, Laptm4a, Atf3
## Gstm2, Csrp1, Ptn, Pltp, Fosb, Dusp1, Igfbp2, Klf4, Krt14, Gas6
## Negative: Tmsb10, Cystm1, Atf5, Plaur, Adam8, Krt6a, Map1b, Tuba1a, Tubb2b, Cst6
## Ptma, Ecm1, Stmn1, Emp3, Uchl1, Upk3bl, Ceacam1, Clcf1, Tm4sf1, Gprc5a
## Cldn7, Tubb6, Lypd3, Cldn3, Lamc2, Phlda2, Ncam1, Cldn4, Mal, Sult2b1
## PC_ 2
## Positive: Igfbp2, Cald1, Tmsb10, Ogn, Tmsb4x, Igfbp4, Adam8, Krt5, Apoe, Plaur
## Fam101b, Tspan6, Abca4, Sparc, Tm4sf1, Col23a1, Lamc2, Krt6a, Clcf1, Tmem108
## Slc29a1, Krt14, Tubb6, Ngf, Cst6, Oaf, Tubb2b, Ecm1, Hmga2, Car12
## Negative: Prdx6, Adh7, Gsta3, Sftpd, Cyp2f2, Ly6a, Tst, Slc16a11, Ifitm1, Reg3g
## Muc1, Por, Cbr2, Adh1, Lrrc26, Selenbp1, Mgst1, Gpx2, Sult1d1, Gstp1
## Cbr3, Cd14, Cyp4b1, Cxcl17, Gsto1, Serpinb11, Elf3, Ly6c1, Serpina3n, Epas1
## PC_ 3
## Positive: Adh7, Aldh3a1, Selenbp1, Adh1, Ly6a, Oat, Defb1, Fmo2, Cbr3, Gstp1
## Tmem176a, Atp1b1, Tmem176b, Efemp1, Igfbp4, Upk1b, Sparc, Ly6c1, Tgm2, Igfbp2
## Mgp, Kcnj16, Sult1d1, Gmpr, Cdo1, Krt15, Nrtn, Cyp4b1, Muc1, Bsg
## Negative: Cldn4, X1810011O10Rik, Klf4, Sfn, Zfp36, S100a10, Mast4, Atf3, Tagln2, Fosb
## Avpi1, Arc, Nr4a1, Cldn7, Nfat5, Egr1, Klf6, Cdkn1a, Phlda1, Lmna
## Crym, Neat1, Anxa1, Krt18, S100a6, Acsm4, Lypd3, Papss2, Ier2, Krt23
## PC_ 4
## Positive: mt.Co1, mt.Co3, mt.Atp6, Hpgd, Pla2r1, X1500009L16Rik, Acsm4, Rps27rt, Car12, Dapl1
## Myo6, Fth1, Papss2, Lypd2, Insl6, Cyb5r3, Scd2, Gas6, Trp63, Nmb
## Eya1, Piezo2, Limd2, Crym, Sema5a, Serpinb1b, Ftl1, X4631405K08Rik, Fndc1, Scgb1c1
## Negative: Actg1, Anxa3, Ier3, Nfkbia, F3, Tmem176b, Krt7, Sfn, Cyr61, Cldn4
## Krt14, Tacstd2, Tuba1c, Tmem176a, Fos, Bhlhe40, Junb, Tnfrsf12a, Gmpr, Cish
## Fst, Cxcl16, Plaur, Krt17, Ubc, Trib1, Oat, Atf3, Adam8, Dusp5
## PC_ 5
## Positive: Fosb, Egr1, Neat1, Irs2, Pdcd4, Nr4a1, Lmna, Klf9, Cbx3, Dst
## Ccdc3, Mast4, Sdc4, Timp3, Nfia, Kdm6b, Birc5, Maff, Spry2, Il6st
## Stmn1, Hist1h2ap, Nfat5, Sik1, Arid5a, Gm26532, Pbk, Gclc, Tpm1, Fam101b
## Negative: Aldh2, Gstm1, Prdx1, Aqp3, Ces1d, Tst, Gsn, Clu, Ephx1, Cpe
## Ifitm3, Ubc, Krt18, Prdx6, Hspb1, Galm, Fmo6, Hspa5, Ppa1, Krt5
## Krt8, Pltp, Mgst1, Ctsl, Cyp2a5, Alas1, Cryl1, Calml4, Gstm2, Actg1rna <- RunUMAP(rna, dims = 1:20, min_dist=.1)## Warning: The following arguments are not used: min_dist## 17:45:23 UMAP embedding parameters a = 0.9922 b = 1.112
## 17:45:23 Read 3064 rows and found 20 numeric columns
## 17:45:23 Using Annoy for neighbor search, n_neighbors = 30
## 17:45:23 Building Annoy index with metric = cosine, n_trees = 50
## 0% 10 20 30 40 50 60 70 80 90 100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 17:45:23 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpBeaM9s/file175b84af1391e
## 17:45:23 Searching Annoy index using 1 thread, search_k = 3000
## 17:45:24 Annoy recall = 100%
## 17:45:25 Commencing smooth kNN distance calibration using 1 thread
## 17:45:27 Initializing from normalized Laplacian + noise
## 17:45:27 Commencing optimization for 500 epochs, with 122146 positive edges
## 17:45:32 Optimization finishedp1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=TRUE) +
NoLegend() +
ggtitle("scATAC-seq")## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.p2 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend() + ggtitle("scRNA-seq")
p1 + p2
# Markers for activated/resting HBC in scRNA-seq
FeaturePlot(rna, features = "Krtdap")
FeaturePlot(rna, features = "Krt6a")
FeaturePlot(rna, features = "Krt5")
FeaturePlot(rna, features = "Krt14")
FeaturePlot(rna, features = "Trp63")
# Markers for activated/resting HBC in scATAC-seq
FeaturePlot(atac, features = "Krtdap") + xlim(c(-5,7))## Warning: Could not find Krtdap in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 14 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Krt6a") + xlim(c(-5,7))## Warning: Could not find Krt6a in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 14 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Krt5") + xlim(c(-6,6))## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Krt14") + xlim(c(-6,6))## Warning: Could not find Krt14 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
FeaturePlot(atac, features = "Trp63") + xlim(c(-6,6))## Warning: Could not find Trp63 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
plots for paper
# RNA
pRNA1 <- DimPlot(rna, group.by = "celltype", label = FALSE, repel = TRUE) +
ggtitle("scRNA-seq: 24h HBC* and rHBC cells") +
scale_color_manual(labels = c("Activated", "Regenerated\n(resting)"),
values = c("steelblue", "salmon")) +
labs(color = "Cell state")
# ggsave("../plots/scRNA24hHBC.png", width=9)
# ATAC
pATAC1 <- DimPlot(atac, reduction = "umap", group.by="clHBCMerged", label=FALSE) +
ggtitle("scATAC-seq HBC cell states") +
xlim(c(-6, 7))
pATAC2 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=FALSE) +
scale_color_manual(labels = c("Injured", "Uninjured"),
values = c("orange", "darkseagreen3")) +
ggtitle("scATAC-seq HBC treatment") +
xlim(c(-6, 7))
cowplot::plot_grid(pATAC1, pATAC2)## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
# ggsave("../plots/ATACClusters_seurat.png", width=10)Data integration
# find anchor features
transfer.anchors <- FindTransferAnchors(reference = rna, query = atac,
features = VariableFeatures(object = rna),
reference.assay = "RNA",
query.assay = "ACTIVITY",
reduction = "cca")## Warning in RunCCA.Seurat(object1 = reference, object2 = query, features =
## features, : Running CCA on different assays## Running CCA## Merging objects## Finding neighborhoods## Finding anchors## Found 10751 anchors## Filtering anchors## Retained 1185 anchors## Extracting within-dataset neighbors# predict ATAC cell type based on RNA labels
celltype.predictions <- TransferData(anchorset = transfer.anchors,
refdata = rna$celltype,
weight.reduction = atac[["lsi"]])## Finding integration vectors## Finding integration vector weights## Predicting cell labels# add predicted labels to ATAC
atac <- AddMetaData(atac, metadata = celltype.predictions)
hist(atac$prediction.score.max)
atac$predicted.id <- factor(atac$predicted.id) # to make the colors match
pat1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=TRUE) +
NoLegend() +
ggtitle("scATAC-seq truth") +
xlim(c(-6,7))
pat2 <- DimPlot(atac, group.by = "predicted.id", label = TRUE, repel = TRUE) +
ggtitle("scATAC-seq predicted") +
NoLegend() +
scale_colour_hue(drop = FALSE) +
xlim(c(-6,7))
pat3 <- DimPlot(atac, group.by = "clHBCMerged", label = TRUE, repel = TRUE) +
ggtitle("scATAC-seq ArchR clusters") +
NoLegend() +
scale_colour_hue(drop = FALSE) +
xlim(c(-6,7))
prna4 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + ggtitle("scRNA-seq") +
NoLegend()
pat1 + pat2 + pat3 + prna4## Warning: Removed 14 rows containing missing values (geom_point).## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
# prediction in ArchR DR space
archrUmap$predicted.id <- celltype.predictions$predicted.id
pat1_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=treatment)) +
geom_point(size = .75) +
ggtitle("scATAC-seq treatment") +
theme_classic() +
scale_color_manual(values = c("orange", "darkseagreen3"))
pat2_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
geom_point(size = .75) +
ggtitle("scATAC-seq predicted") +
theme_classic()
pat3_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=clHBCMerged)) +
geom_point(size = .75) +
ggtitle("scATAC-seq clusters") +
theme_classic()
prna1_2 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) +
ggtitle("scRNA-seq") +
NoLegend()
pat1_2 + pat2_2 + pat3_2 + prna1_2
Compare predictions with truth and ArchR clusters
library(ggalluvial)
dfCellWide <- data.frame(treatment=atac$treatment,
predicted=atac$predicted.id,
cluster= archProj@cellColData$clHBCMerged)
pAlluv <- ggplot(dfCellWide,
aes( axis1 = cluster, axis2 = predicted, axis3=treatment)) +
geom_alluvium(aes(fill=treatment), width = 1/12) +
scale_x_discrete(limits = c("Cluster", "Predicted", "Treatment"), expand = c(.05, .05)) +
geom_stratum(width = 1/12, fill = "white", color = "grey") +
geom_label(stat = "stratum", aes(label = after_stat(stratum))) +
theme(panel.background = element_rect(fill="white"))
pAlluv
Check gene expression of chromatin-responsive genes upon injury
Here, we check the gene expression for genes that were found to be DA in the scATAC-seq data when testing for the injury effect for both the injured vs uninjured cluster and injury vs uninjured cells in the hybrid cluster.
Chromatin more accessible in injury
chromUpInjury <- read.table("../data/scATAC/sharedUpInjury.txt",
stringsAsFactors = FALSE)[,1]
chromUpInjury <- chromUpInjury[chromUpInjury %in% rownames(countHBC)]
plotByExpr <- function(rna, genes, counts, name="feature", scale=TRUE){
genes <- genes[genes %in% rownames(countHBC)]
if(scale){
scaledCounts <- t(scale(t(counts[genes,]+1e-5)))
scaledCounts <- scaledCounts[!is.na(rowVars(scaledCounts)),]
yhat <- colSums(scaledCounts)
} else {
yhat <- colSums(counts[genes,]+1e-5)
}
rna <- AddMetaData(rna, yhat, col.name = name)
FeaturePlot(rna, features = name) + theme(legend.position = "none")
}
plotByExpr(rna, chromUpInjury, countHBC, "chromUpInjurySum", scale=FALSE)
plotByExpr(rna, chromUpInjury, countHBC, "chromUpInjurySum_scaled", scale=TRUE)
Chromatin less accessible in injury
chromDownInjury <- read.table("../data/scATAC/sharedDownInjury.txt",
stringsAsFactors = FALSE)[,1]
chromDownInjury <- chromDownInjury[chromDownInjury %in% rownames(countHBC)]
plotByExpr(rna, chromDownInjury, countHBC, "chromDownInjurySum", scale=FALSE)
plotByExpr(rna, chromDownInjury, countHBC, "chromDownjurySum_scaled", scale=TRUE)
Check gene expression of chromatin-responsive genes in hybrid vs uninjured ATAC-seq clusters
chromHybrid <- read.table("../data/scATAC/sharedHybridGenes.txt",
stringsAsFactors = FALSE)[,1]
chromHybrid <- chromHybrid[chromHybrid %in% rownames(countHBC)]
plotByExpr(rna, chromHybrid, countHBC, "hybridSum", scale=FALSE)
plotByExpr(rna, chromHybrid, countHBC, "chybridSum_scaled", scale=TRUE)
ATAC-seq cluster markers
# atacMarkers <- readRDS("../../../scATAC/data/markerList_atac_ArchR_hbcClusters.rds")
atacMarkers <- readRDS("../data/scATAC/markerList_atac_ArchR_hbcClusters.rds")
pc1 <- plotByExpr(rna, atacMarkers$C1_Inj$name, countHBC, "C1_Inj", scale=TRUE)
pc2 <- plotByExpr(rna, atacMarkers$C2_Hybrid$name, countHBC, "C2_Hybrid", scale=TRUE)
pc3 <- plotByExpr(rna, atacMarkers$C3_UI$name, countHBC, "C3_UI", scale=TRUE)
pMarkersATAC <- cowplot::plot_grid(plotlist=list(pc1,pc2,pc3), nrow=1)
# ggsave("../plots/ATACMarkersOnRNA.png", width=19)
## are the ATAC markers driven by both injured and uninjured cells in hybrid cluster? Yes.
activityHybridMarkers <- as.matrix(atac[["ACTIVITY"]]@counts)[atacMarkers$C2_Hybrid$name,]
df <- data.frame(y=colSums(t(scale(t(activityHybridMarkers)))),
state=factor(interaction(atac$clHBCMerged, atac$treatment)))
ggplot(df, aes(x=state, y=y)) +
geom_boxplot() +
xlab("Group (cell state x treatment)") +
ylab("Scaled activity of hybrid markers") +
ggtitle("Cell state expression of hybrid markers") +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5))
# ggsave("../plots/atacHybridMarkerCellContribution.png", width=9)Redo cell label prediction when assigning hybrid labels
For reference, old plot
# Seurat
pat1 + pat2 + pat3 + prna4## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
# ArchR DR space
pat1_2 + pat2_2 + pat3_2 + prna1_2
Per hybrid cluster separately
ctHybridSub <- as.character(rna$celltype)
hlpid <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1)
ctHybridSub[hlpid] <- NA
hlp2id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1 & rna@reductions$umap@cell.embeddings[,2]<4)
ctHybridSub[hlp2id] <- "hybrid1"
hlp1id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>4)
ctHybridSub[hlp1id] <- "hybrid2"
rna$ctHybridSub <- ctHybridSub
DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE) +
ggtitle("scRNA-seq") 
# predict ATAC cell type based on RNA labels
celltype.predictions3 <- TransferData(anchorset = transfer.anchors,
refdata = rna$ctHybridSub,
weight.reduction = atac[["lsi"]])## Finding integration vectors## Finding integration vector weights## Predicting cell labels# add predicted labels to ATAC
atac$predictedCtHybrid <- celltype.predictions3$predicted.id
atac$predScore <- celltype.predictions3$prediction.score.max
atac$predicted.id <- factor(atac$predicted.id) # to make the colors match
# prediction in Seurat RD space
ppred1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=FALSE) +
ggtitle("scATAC-seq treatment") +
xlim(c(-6, 7))
ppred2 <- DimPlot(atac, group.by = "predictedCtHybrid", label = FALSE, repel = TRUE) +
ggtitle("scATAC-seq predicted state") +
scale_colour_hue(drop = FALSE) +
xlim(c(-6, 7))
ppred3 <- DimPlot(atac, group.by = "clHBCMerged", label = FALSE, repel = TRUE) +
ggtitle("scATAC-seq cell state") +
scale_colour_hue(drop = FALSE) +
xlim(c(-6, 7))
ppred4 <- DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE) +
ggtitle("scRNA-seq")
ppred1 + ppred2 + ppred3 + ppred4## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
## Warning: Removed 14 rows containing missing values (geom_point).
# ggsave("../plots/cellTransferSubHybrid.png", width=10)
## size of cell is prediction prob
DimPlot(atac, group.by = "predictedCtHybrid", label = TRUE, repel = TRUE, pt.size=celltype.predictions3$prediction.score.max*1.25) + ggtitle("scATAC-seq predicted") +
NoLegend() + scale_colour_hue(drop = FALSE) + xlim(c(-6,6))## Warning: Removed 55 rows containing missing values (geom_point).
boxplot(celltype.predictions3$prediction.score.max ~ celltype.predictions3$predicted.id,
xlab="Predicted state", ylab="Probability")
# prediction in ArchR DR space
archrUmap$predicted.id <- celltype.predictions3$predicted.id
ppred1_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=treatment)) +
geom_point(size = .75) +
ggtitle("scATAC-seq treatment") +
theme_classic() +
scale_color_manual(values = c("orange", "darkseagreen3"))
ppred2_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
geom_point(size = .5) +
ggtitle("scATAC-seq predicted") +
theme_classic()
ppred3_2 <- ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=clHBCMerged)) +
geom_point(size = .75) +
ggtitle("scATAC-seq clusters") +
theme_classic()
ppred4_2 <- DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE, pt.size=.1) +
ggtitle("scRNA-seq") +
theme(axis.title = element_text(size = 11),
axis.text = element_text(size = 9),
plot.title = element_text(face = "plain", size=13))
ppred1_2 + ppred2_2 + ppred3_2 + ppred4_2
## size of cell is prediction prob
ggplot(archrUmap, aes(x=UMAP1, y=UMAP2, color=predicted.id)) +
geom_point(size = celltype.predictions3$prediction.score.max*1.25) +
ggtitle("scATAC-seq predicted") +
theme_classic()
boxplot(celltype.predictions3$prediction.score.max ~ celltype.predictions3$predicted.id,
xlab="Predicted state", ylab="Probability")
Markers for sub-hybrid clusters
## add annotation according to prediction
predictedAnn <- archProj@cellColData$clHBCMerged
predictedAnn[predictedAnn == "C3_UI"] <- "resting"
predictedAnn[predictedAnn == "C1_inj"] <- "activated"
hybId <- predictedAnn == "C2_hybrid"
predictedAnn[which(hybId)] <- NA
predictedAnn[hybId &
!celltype.predictions3$predicted.id == "hybrid1" &
!celltype.predictions3$predicted.id == "hybrid2"] <- "hybrid"
predictedAnn[which(celltype.predictions3$predicted.id == "hybrid1")] <- "hybrid1"
predictedAnn[which(celltype.predictions3$predicted.id == "hybrid2")] <- "hybrid2"
table(predictedAnn)## predictedAnn
## C1_Inj C2_Hybrid hybrid1 hybrid2 resting
## 1485 933 197 223 1238## visualize with Seurat
atac$predAnn <- predictedAnn
DimPlot(atac, group.by = "predAnn", label = TRUE, repel = TRUE) + ggtitle("scATAC-seq predicted") +
NoLegend() + scale_colour_hue(drop = FALSE) + xlim(c(-6,7))## Warning: Removed 14 rows containing missing values (geom_point).
## distribution of injured and uninjured cells in each cluster
df <- data.frame(treatment=atac$treatment,
cluster=atac$predAnn)
dfTotal <- df %>%
group_by(cluster) %>%
mutate(totalCells=n()) %>%
group_by(cluster, treatment) %>%
summarize(fracCells = n() / unique(totalCells))
ggplot(dfTotal, aes(x=treatment, fill=treatment, y=fracCells)) +
geom_bar(position="dodge", stat="identity") +
facet_wrap(.~cluster) +
ylab("Fraction of cells")
# saveRDS(celltype.predictions3, file="../data/scATAC/20201105_celltypePredictions3.rds")
## are the ATAC markers driven by both injured and uninjured cells in all hybrid cluster? Yes.
activityHybridMarkers <- as.matrix(atac[["ACTIVITY"]]@counts)[atacMarkers$C2_Hybrid$name,]
df <- data.frame(y=colSums(t(scale(t(activityHybridMarkers)))),
state=factor(interaction(atac$predAnn, atac$treatment)))
ggplot(df, aes(x=state, y=y)) +
geom_boxplot() +
xlab("Group (cell state x treatment)") +
ylab("Scaled activity of hybrid markers") +
ggtitle("Cell state expression of hybrid markers") +
theme_bw() +
theme(axis.text.x = element_text(angle=45, vjust=0.5))
# ggsave("../plots/atacHybridMarkerCellContributionAllHybridClusters.png", width=9)New peak calling based on sub-clustering of hybrid cells as suggested by RNA
archProj@cellColData$predAnn <- predictedAnn
# pseudobulk for peak calling
archProj <- addGroupCoverages(ArchRProj = archProj, groupBy = "predAnn",
minReplicates = 2, maxReplicates = 6,
force = TRUE,
verbose = FALSE)## ArchR logging to : ArchRLogs/ArchR-addGroupCoverages-175b878f01d6b-Date-2020-12-03_Time-17-47-32.log
## If there is an issue, please report to github with logFile!## C1_Inj (1 of 5) : CellGroups N = 4## C2_Hybrid (2 of 5) : CellGroups N = 6## hybrid1 (3 of 5) : CellGroups N = 3## hybrid2 (4 of 5) : CellGroups N = 2## resting (5 of 5) : CellGroups N = 5## 2020-12-03 17:47:33 : Creating Coverage Files!, 0.027 mins elapsed.## 2020-12-03 17:47:33 : Batch Execution w/ safelapply!, 0.027 mins elapsed.## Number of Cells = 500## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 263## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 163## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 54## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 269## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 252## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 140## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 113## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 83## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 76## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 94## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 57## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 46## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 170## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 53## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 500## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 355## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 178## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 54## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## Number of Cells = 51## Coverage File Exists!## Added Coverage Group## Added Metadata Group## Added ArrowCoverage Class## Added Coverage/Info## Added Coverage/Info/CellNames## 2020-12-03 17:55:07 : Adding Kmer Bias to Coverage Files!, 7.589 mins elapsed.## Kmer Bias chr1 (1 of 21)## chr1 Coverage File chr1 (1 of 20)
## Coverage File chr1 (2 of 20)
## Coverage File chr1 (3 of 20)
## Coverage File chr1 (4 of 20)
## Coverage File chr1 (5 of 20)
## Coverage File chr1 (6 of 20)
## Coverage File chr1 (7 of 20)
## Coverage File chr1 (8 of 20)
## Coverage File chr1 (9 of 20)
## Coverage File chr1 (10 of 20)
## Coverage File chr1 (11 of 20)
## Coverage File chr1 (12 of 20)
## Coverage File chr1 (13 of 20)
## Coverage File chr1 (14 of 20)
## Coverage File chr1 (15 of 20)
## Coverage File chr1 (16 of 20)
## Coverage File chr1 (17 of 20)
## Coverage File chr1 (18 of 20)
## Coverage File chr1 (19 of 20)
## Coverage File chr1 (20 of 20)
## Kmer Bias chr10 (2 of 21)
## chr10 Coverage File chr10 (1 of 20)
## Coverage File chr10 (2 of 20)
## Coverage File chr10 (3 of 20)
## Coverage File chr10 (4 of 20)
## Coverage File chr10 (5 of 20)
## Coverage File chr10 (6 of 20)
## Coverage File chr10 (7 of 20)
## Coverage File chr10 (8 of 20)
## Coverage File chr10 (9 of 20)
## Coverage File chr10 (10 of 20)
## Coverage File chr10 (11 of 20)
## Coverage File chr10 (12 of 20)
## Coverage File chr10 (13 of 20)
## Coverage File chr10 (14 of 20)
## Coverage File chr10 (15 of 20)
## Coverage File chr10 (16 of 20)
## Coverage File chr10 (17 of 20)
## Coverage File chr10 (18 of 20)
## Coverage File chr10 (19 of 20)
## Coverage File chr10 (20 of 20)
## Kmer Bias chr11 (3 of 21)
## chr11 Coverage File chr11 (1 of 20)
## Coverage File chr11 (2 of 20)
## Coverage File chr11 (3 of 20)
## Coverage File chr11 (4 of 20)
## Coverage File chr11 (5 of 20)
## Coverage File chr11 (6 of 20)
## Coverage File chr11 (7 of 20)
## Coverage File chr11 (8 of 20)
## Coverage File chr11 (9 of 20)
## Coverage File chr11 (10 of 20)
## Coverage File chr11 (11 of 20)
## Coverage File chr11 (12 of 20)
## Coverage File chr11 (13 of 20)
## Coverage File chr11 (14 of 20)
## Coverage File chr11 (15 of 20)
## Coverage File chr11 (16 of 20)
## Coverage File chr11 (17 of 20)
## Coverage File chr11 (18 of 20)
## Coverage File chr11 (19 of 20)
## Coverage File chr11 (20 of 20)
## Kmer Bias chr12 (4 of 21)
## chr12 Coverage File chr12 (1 of 20)
## Coverage File chr12 (2 of 20)
## Coverage File chr12 (3 of 20)
## Coverage File chr12 (4 of 20)
## Coverage File chr12 (5 of 20)
## Coverage File chr12 (6 of 20)
## Coverage File chr12 (7 of 20)
## Coverage File chr12 (8 of 20)
## Coverage File chr12 (9 of 20)
## Coverage File chr12 (10 of 20)
## Coverage File chr12 (11 of 20)
## Coverage File chr12 (12 of 20)
## Coverage File chr12 (13 of 20)
## Coverage File chr12 (14 of 20)
## Coverage File chr12 (15 of 20)
## Coverage File chr12 (16 of 20)
## Coverage File chr12 (17 of 20)
## Coverage File chr12 (18 of 20)
## Coverage File chr12 (19 of 20)
## Coverage File chr12 (20 of 20)
## Kmer Bias chr13 (5 of 21)
## chr13 Coverage File chr13 (1 of 20)
## Coverage File chr13 (2 of 20)
## Coverage File chr13 (3 of 20)
## Coverage File chr13 (4 of 20)
## Coverage File chr13 (5 of 20)
## Coverage File chr13 (6 of 20)
## Coverage File chr13 (7 of 20)
## Coverage File chr13 (8 of 20)
## Coverage File chr13 (9 of 20)
## Coverage File chr13 (10 of 20)
## Coverage File chr13 (11 of 20)
## Coverage File chr13 (12 of 20)
## Coverage File chr13 (13 of 20)
## Coverage File chr13 (14 of 20)
## Coverage File chr13 (15 of 20)
## Coverage File chr13 (16 of 20)
## Coverage File chr13 (17 of 20)
## Coverage File chr13 (18 of 20)
## Coverage File chr13 (19 of 20)
## Coverage File chr13 (20 of 20)
## Kmer Bias chr14 (6 of 21)
## chr14 Coverage File chr14 (1 of 20)
## Coverage File chr14 (2 of 20)
## Coverage File chr14 (3 of 20)
## Coverage File chr14 (4 of 20)
## Coverage File chr14 (5 of 20)
## Coverage File chr14 (6 of 20)
## Coverage File chr14 (7 of 20)
## Coverage File chr14 (8 of 20)
## Coverage File chr14 (9 of 20)
## Coverage File chr14 (10 of 20)
## Coverage File chr14 (11 of 20)
## Coverage File chr14 (12 of 20)
## Coverage File chr14 (13 of 20)
## Coverage File chr14 (14 of 20)
## Coverage File chr14 (15 of 20)
## Coverage File chr14 (16 of 20)
## Coverage File chr14 (17 of 20)
## Coverage File chr14 (18 of 20)
## Coverage File chr14 (19 of 20)
## Coverage File chr14 (20 of 20)
## Kmer Bias chr15 (7 of 21)
## chr15 Coverage File chr15 (1 of 20)
## Coverage File chr15 (2 of 20)
## Coverage File chr15 (3 of 20)
## Coverage File chr15 (4 of 20)
## Coverage File chr15 (5 of 20)
## Coverage File chr15 (6 of 20)
## Coverage File chr15 (7 of 20)
## Coverage File chr15 (8 of 20)
## Coverage File chr15 (9 of 20)
## Coverage File chr15 (10 of 20)
## Coverage File chr15 (11 of 20)
## Coverage File chr15 (12 of 20)
## Coverage File chr15 (13 of 20)
## Coverage File chr15 (14 of 20)
## Coverage File chr15 (15 of 20)
## Coverage File chr15 (16 of 20)
## Coverage File chr15 (17 of 20)
## Coverage File chr15 (18 of 20)
## Coverage File chr15 (19 of 20)
## Coverage File chr15 (20 of 20)
## Kmer Bias chr16 (8 of 21)
## chr16 Coverage File chr16 (1 of 20)
## Coverage File chr16 (2 of 20)
## Coverage File chr16 (3 of 20)
## Coverage File chr16 (4 of 20)
## Coverage File chr16 (5 of 20)
## Coverage File chr16 (6 of 20)
## Coverage File chr16 (7 of 20)
## Coverage File chr16 (8 of 20)
## Coverage File chr16 (9 of 20)
## Coverage File chr16 (10 of 20)
## Coverage File chr16 (11 of 20)
## Coverage File chr16 (12 of 20)
## Coverage File chr16 (13 of 20)
## Coverage File chr16 (14 of 20)
## Coverage File chr16 (15 of 20)
## Coverage File chr16 (16 of 20)
## Coverage File chr16 (17 of 20)
## Coverage File chr16 (18 of 20)
## Coverage File chr16 (19 of 20)
## Coverage File chr16 (20 of 20)
## Kmer Bias chr17 (9 of 21)
## chr17 Coverage File chr17 (1 of 20)
## Coverage File chr17 (2 of 20)
## Coverage File chr17 (3 of 20)
## Coverage File chr17 (4 of 20)
## Coverage File chr17 (5 of 20)
## Coverage File chr17 (6 of 20)
## Coverage File chr17 (7 of 20)
## Coverage File chr17 (8 of 20)
## Coverage File chr17 (9 of 20)
## Coverage File chr17 (10 of 20)
## Coverage File chr17 (11 of 20)
## Coverage File chr17 (12 of 20)
## Coverage File chr17 (13 of 20)
## Coverage File chr17 (14 of 20)
## Coverage File chr17 (15 of 20)
## Coverage File chr17 (16 of 20)
## Coverage File chr17 (17 of 20)
## Coverage File chr17 (18 of 20)
## Coverage File chr17 (19 of 20)
## Coverage File chr17 (20 of 20)
## Kmer Bias chr18 (10 of 21)
## chr18 Coverage File chr18 (1 of 20)
## Coverage File chr18 (2 of 20)
## Coverage File chr18 (3 of 20)
## Coverage File chr18 (4 of 20)
## Coverage File chr18 (5 of 20)
## Coverage File chr18 (6 of 20)
## Coverage File chr18 (7 of 20)
## Coverage File chr18 (8 of 20)
## Coverage File chr18 (9 of 20)
## Coverage File chr18 (10 of 20)
## Coverage File chr18 (11 of 20)
## Coverage File chr18 (12 of 20)
## Coverage File chr18 (13 of 20)
## Coverage File chr18 (14 of 20)
## Coverage File chr18 (15 of 20)
## Coverage File chr18 (16 of 20)
## Coverage File chr18 (17 of 20)
## Coverage File chr18 (18 of 20)
## Coverage File chr18 (19 of 20)
## Coverage File chr18 (20 of 20)
## Kmer Bias chr19 (11 of 21)
## chr19 Coverage File chr19 (1 of 20)
## Coverage File chr19 (2 of 20)
## Coverage File chr19 (3 of 20)
## Coverage File chr19 (4 of 20)
## Coverage File chr19 (5 of 20)
## Coverage File chr19 (6 of 20)
## Coverage File chr19 (7 of 20)
## Coverage File chr19 (8 of 20)
## Coverage File chr19 (9 of 20)
## Coverage File chr19 (10 of 20)
## Coverage File chr19 (11 of 20)
## Coverage File chr19 (12 of 20)
## Coverage File chr19 (13 of 20)
## Coverage File chr19 (14 of 20)
## Coverage File chr19 (15 of 20)
## Coverage File chr19 (16 of 20)
## Coverage File chr19 (17 of 20)
## Coverage File chr19 (18 of 20)
## Coverage File chr19 (19 of 20)
## Coverage File chr19 (20 of 20)
## Kmer Bias chr2 (12 of 21)
## chr2 Coverage File chr2 (1 of 20)
## Coverage File chr2 (2 of 20)
## Coverage File chr2 (3 of 20)
## Coverage File chr2 (4 of 20)
## Coverage File chr2 (5 of 20)
## Coverage File chr2 (6 of 20)
## Coverage File chr2 (7 of 20)
## Coverage File chr2 (8 of 20)
## Coverage File chr2 (9 of 20)
## Coverage File chr2 (10 of 20)
## Coverage File chr2 (11 of 20)
## Coverage File chr2 (12 of 20)
## Coverage File chr2 (13 of 20)
## Coverage File chr2 (14 of 20)
## Coverage File chr2 (15 of 20)
## Coverage File chr2 (16 of 20)
## Coverage File chr2 (17 of 20)
## Coverage File chr2 (18 of 20)
## Coverage File chr2 (19 of 20)
## Coverage File chr2 (20 of 20)
## Kmer Bias chr3 (13 of 21)
## chr3 Coverage File chr3 (1 of 20)
## Coverage File chr3 (2 of 20)
## Coverage File chr3 (3 of 20)
## Coverage File chr3 (4 of 20)
## Coverage File chr3 (5 of 20)
## Coverage File chr3 (6 of 20)
## Coverage File chr3 (7 of 20)
## Coverage File chr3 (8 of 20)
## Coverage File chr3 (9 of 20)
## Coverage File chr3 (10 of 20)
## Coverage File chr3 (11 of 20)
## Coverage File chr3 (12 of 20)
## Coverage File chr3 (13 of 20)
## Coverage File chr3 (14 of 20)
## Coverage File chr3 (15 of 20)
## Coverage File chr3 (16 of 20)
## Coverage File chr3 (17 of 20)
## Coverage File chr3 (18 of 20)
## Coverage File chr3 (19 of 20)
## Coverage File chr3 (20 of 20)
## Kmer Bias chr4 (14 of 21)
## chr4 Coverage File chr4 (1 of 20)
## Coverage File chr4 (2 of 20)
## Coverage File chr4 (3 of 20)
## Coverage File chr4 (4 of 20)
## Coverage File chr4 (5 of 20)
## Coverage File chr4 (6 of 20)
## Coverage File chr4 (7 of 20)
## Coverage File chr4 (8 of 20)
## Coverage File chr4 (9 of 20)
## Coverage File chr4 (10 of 20)
## Coverage File chr4 (11 of 20)
## Coverage File chr4 (12 of 20)
## Coverage File chr4 (13 of 20)
## Coverage File chr4 (14 of 20)
## Coverage File chr4 (15 of 20)
## Coverage File chr4 (16 of 20)
## Coverage File chr4 (17 of 20)
## Coverage File chr4 (18 of 20)
## Coverage File chr4 (19 of 20)
## Coverage File chr4 (20 of 20)
## Kmer Bias chr5 (15 of 21)
## chr5 Coverage File chr5 (1 of 20)
## Coverage File chr5 (2 of 20)
## Coverage File chr5 (3 of 20)
## Coverage File chr5 (4 of 20)
## Coverage File chr5 (5 of 20)
## Coverage File chr5 (6 of 20)
## Coverage File chr5 (7 of 20)
## Coverage File chr5 (8 of 20)
## Coverage File chr5 (9 of 20)
## Coverage File chr5 (10 of 20)
## Coverage File chr5 (11 of 20)
## Coverage File chr5 (12 of 20)
## Coverage File chr5 (13 of 20)
## Coverage File chr5 (14 of 20)
## Coverage File chr5 (15 of 20)
## Coverage File chr5 (16 of 20)
## Coverage File chr5 (17 of 20)
## Coverage File chr5 (18 of 20)
## Coverage File chr5 (19 of 20)
## Coverage File chr5 (20 of 20)
## Kmer Bias chr6 (16 of 21)
## chr6 Coverage File chr6 (1 of 20)
## Coverage File chr6 (2 of 20)
## Coverage File chr6 (3 of 20)
## Coverage File chr6 (4 of 20)
## Coverage File chr6 (5 of 20)
## Coverage File chr6 (6 of 20)
## Coverage File chr6 (7 of 20)
## Coverage File chr6 (8 of 20)
## Coverage File chr6 (9 of 20)
## Coverage File chr6 (10 of 20)
## Coverage File chr6 (11 of 20)
## Coverage File chr6 (12 of 20)
## Coverage File chr6 (13 of 20)
## Coverage File chr6 (14 of 20)
## Coverage File chr6 (15 of 20)
## Coverage File chr6 (16 of 20)
## Coverage File chr6 (17 of 20)
## Coverage File chr6 (18 of 20)
## Coverage File chr6 (19 of 20)
## Coverage File chr6 (20 of 20)
## Kmer Bias chr7 (17 of 21)
## chr7 Coverage File chr7 (1 of 20)
## Coverage File chr7 (2 of 20)
## Coverage File chr7 (3 of 20)
## Coverage File chr7 (4 of 20)
## Coverage File chr7 (5 of 20)
## Coverage File chr7 (6 of 20)
## Coverage File chr7 (7 of 20)
## Coverage File chr7 (8 of 20)
## Coverage File chr7 (9 of 20)
## Coverage File chr7 (10 of 20)
## Coverage File chr7 (11 of 20)
## Coverage File chr7 (12 of 20)
## Coverage File chr7 (13 of 20)
## Coverage File chr7 (14 of 20)
## Coverage File chr7 (15 of 20)
## Coverage File chr7 (16 of 20)
## Coverage File chr7 (17 of 20)
## Coverage File chr7 (18 of 20)
## Coverage File chr7 (19 of 20)
## Coverage File chr7 (20 of 20)
## Kmer Bias chr8 (18 of 21)
## chr8 Coverage File chr8 (1 of 20)
## Coverage File chr8 (2 of 20)
## Coverage File chr8 (3 of 20)
## Coverage File chr8 (4 of 20)
## Coverage File chr8 (5 of 20)
## Coverage File chr8 (6 of 20)
## Coverage File chr8 (7 of 20)
## Coverage File chr8 (8 of 20)
## Coverage File chr8 (9 of 20)
## Coverage File chr8 (10 of 20)
## Coverage File chr8 (11 of 20)
## Coverage File chr8 (12 of 20)
## Coverage File chr8 (13 of 20)
## Coverage File chr8 (14 of 20)
## Coverage File chr8 (15 of 20)
## Coverage File chr8 (16 of 20)
## Coverage File chr8 (17 of 20)
## Coverage File chr8 (18 of 20)
## Coverage File chr8 (19 of 20)
## Coverage File chr8 (20 of 20)
## Kmer Bias chr9 (19 of 21)
## chr9 Coverage File chr9 (1 of 20)
## Coverage File chr9 (2 of 20)
## Coverage File chr9 (3 of 20)
## Coverage File chr9 (4 of 20)
## Coverage File chr9 (5 of 20)
## Coverage File chr9 (6 of 20)
## Coverage File chr9 (7 of 20)
## Coverage File chr9 (8 of 20)
## Coverage File chr9 (9 of 20)
## Coverage File chr9 (10 of 20)
## Coverage File chr9 (11 of 20)
## Coverage File chr9 (12 of 20)
## Coverage File chr9 (13 of 20)
## Coverage File chr9 (14 of 20)
## Coverage File chr9 (15 of 20)
## Coverage File chr9 (16 of 20)
## Coverage File chr9 (17 of 20)
## Coverage File chr9 (18 of 20)
## Coverage File chr9 (19 of 20)
## Coverage File chr9 (20 of 20)
## Kmer Bias chrX (20 of 21)
## chrX Coverage File chrX (1 of 20)
## Coverage File chrX (2 of 20)
## Coverage File chrX (3 of 20)
## Coverage File chrX (4 of 20)
## Coverage File chrX (5 of 20)
## Coverage File chrX (6 of 20)
## Coverage File chrX (7 of 20)
## Coverage File chrX (8 of 20)
## Coverage File chrX (9 of 20)
## Coverage File chrX (10 of 20)
## Coverage File chrX (11 of 20)
## Coverage File chrX (12 of 20)
## Coverage File chrX (13 of 20)
## Coverage File chrX (14 of 20)
## Coverage File chrX (15 of 20)
## Coverage File chrX (16 of 20)
## Coverage File chrX (17 of 20)
## Coverage File chrX (18 of 20)
## Coverage File chrX (19 of 20)
## Coverage File chrX (20 of 20)
## Kmer Bias chrY (21 of 21)
## chrY Coverage File chrY (1 of 20)
## Coverage File chrY (2 of 20)
## Coverage File chrY (3 of 20)
## Coverage File chrY (4 of 20)
## Coverage File chrY (5 of 20)
## Coverage File chrY (6 of 20)
## Coverage File chrY (7 of 20)
## Coverage File chrY (8 of 20)
## Coverage File chrY (9 of 20)
## Coverage File chrY (10 of 20)
## Coverage File chrY (11 of 20)
## Coverage File chrY (12 of 20)
## Coverage File chrY (13 of 20)
## Coverage File chrY (14 of 20)
## Coverage File chrY (15 of 20)
## Coverage File chrY (16 of 20)
## Coverage File chrY (17 of 20)
## Coverage File chrY (18 of 20)
## Coverage File chrY (19 of 20)
## Coverage File chrY (20 of 20)
## Completed Kmer Bias Calculation
## Adding Kmer Bias (1 of 20)
## Adding Kmer Bias (2 of 20)
## Adding Kmer Bias (3 of 20)
## Adding Kmer Bias (4 of 20)
## Adding Kmer Bias (5 of 20)
## Adding Kmer Bias (6 of 20)
## Adding Kmer Bias (7 of 20)
## Adding Kmer Bias (8 of 20)
## Adding Kmer Bias (9 of 20)
## Adding Kmer Bias (10 of 20)
## Adding Kmer Bias (11 of 20)
## Adding Kmer Bias (12 of 20)
## Adding Kmer Bias (13 of 20)
## Adding Kmer Bias (14 of 20)
## Adding Kmer Bias (15 of 20)
## Adding Kmer Bias (16 of 20)
## Adding Kmer Bias (17 of 20)
## Adding Kmer Bias (18 of 20)
## Adding Kmer Bias (19 of 20)
## Adding Kmer Bias (20 of 20)
## 2020-12-03 18:02:21 : Finished Creation of Coverage Files!, 14.824 mins elapsed.
## ArchR logging successful to : ArchRLogs/ArchR-addGroupCoverages-175b878f01d6b-Date-2020-12-03_Time-17-47-32.log# peak calling
archProj <- addReproduciblePeakSet(
ArchRProj = archProj,
groupBy = "predAnn",
pathToMacs2 = "/Users/koenvandenberge/opt/anaconda3/bin/macs2",
verbose = FALSE)## ArchR logging to : ArchRLogs/ArchR-addReproduciblePeakSet-175b85c1b99d4-Date-2020-12-03_Time-18-02-21.log
## If there is an issue, please report to github with logFile!
## Calling Peaks with Macs2
## 2020-12-03 18:02:22 : Batch Execution w/ safelapply!, 0 mins elapsed.
## 2020-12-03 18:02:22 : Group 1 of 20, Calling Peaks with MACS2!, 0 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s1Inj-1 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s1Inj-1.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:04:08 : Group 2 of 20, Calling Peaks with MACS2!, 1.765 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s2Inj-2 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s2Inj-2.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:05:42 : Group 3 of 20, Calling Peaks with MACS2!, 3.341 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.s3Inj-3 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.s3Inj-3.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:07:06 : Group 4 of 20, Calling Peaks with MACS2!, 4.731 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C1_Inj._.Other-4 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C1_Inj._.Other-4.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:08:00 : Group 5 of 20, Calling Peaks with MACS2!, 5.644 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s1Inj-5 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s1Inj-5.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:09:09 : Group 6 of 20, Calling Peaks with MACS2!, 6.788 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s3Inj-6 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s3Inj-6.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:10:52 : Group 7 of 20, Calling Peaks with MACS2!, 8.507 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s2UI-7 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s2UI-7.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:11:41 : Group 8 of 20, Calling Peaks with MACS2!, 9.322 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s3UI-8 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s3UI-8.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:12:29 : Group 9 of 20, Calling Peaks with MACS2!, 10.114 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s2Inj-9 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s2Inj-9.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:13:22 : Group 10 of 20, Calling Peaks with MACS2!, 11 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name C2_Hybrid._.s1UI-10 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/C2_Hybrid._.s1UI-10.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:14:17 : Group 11 of 20, Calling Peaks with MACS2!, 11.92 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.s3UI-11 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.s3UI-11.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:15:32 : Group 12 of 20, Calling Peaks with MACS2!, 13.165 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.s3Inj-12 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.s3Inj-12.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:16:37 : Group 13 of 20, Calling Peaks with MACS2!, 14.258 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid1._.Other-13 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid1._.Other-13.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:17:42 : Group 14 of 20, Calling Peaks with MACS2!, 15.334 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid2._.s3UI-14 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid2._.s3UI-14.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:19:11 : Group 15 of 20, Calling Peaks with MACS2!, 16.816 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name hybrid2._.1-15 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/hybrid2._.1-15.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:19:47 : Group 16 of 20, Calling Peaks with MACS2!, 17.426 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s2UI-16 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s2UI-16.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:21:20 : Group 17 of 20, Calling Peaks with MACS2!, 18.968 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s1UI-17 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s1UI-17.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:22:37 : Group 18 of 20, Calling Peaks with MACS2!, 20.256 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s3UI-18 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s3UI-18.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:23:39 : Group 19 of 20, Calling Peaks with MACS2!, 21.284 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.s3Inj-19 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.s3Inj-19.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## 2020-12-03 18:24:20 : Group 20 of 20, Calling Peaks with MACS2!, 21.966 mins elapsed.
## Running Macs2 with Params : macs2 callpeak -g 1.87e+09 --name resting._.Other-20 --treatment /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds/resting._.Other-20.insertions.bed --outdir /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/InsertionBeds --format BED --call-summits --keep-dup all --nomodel --nolambda --shift -75 --extsize 150 -q 0.1
## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/C1_Inj-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/C2_Hybrid-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/hybrid1-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/hybrid2-reproduciblePeaks.gr.rds"## Annotating Peaks : Nearest Gene
## Annotating Peaks : Gene
## Annotating Peaks : TSS
## Annotating Peaks : GC## [1] "/Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/data/scATAC/archR_samples_pairwise_v2/PeakCalls/resting-reproduciblePeaks.gr.rds"## Converged after 9 iterations!
## Plotting Ggplot!archProj <- addPeakMatrix(archProj)## ArchR logging to : ArchRLogs/ArchR-addPeakMatrix-175b85a2b745d-Date-2020-12-03_Time-18-25-47.log
## If there is an issue, please report to github with logFile!
## 2020-12-03 18:25:47 : Batch Execution w/ safelapply!, 0 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2020-12-03 18:25:55 : Adding s1Inj to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2020-12-03 18:25:57 : Adding s1Inj to PeakMatrix for Chr (2 of 20)!, 0.054 mins elapsed.
## 2020-12-03 18:26:00 : Adding s1Inj to PeakMatrix for Chr (3 of 20)!, 0.097 mins elapsed.
## 2020-12-03 18:26:02 : Adding s1Inj to PeakMatrix for Chr (4 of 20)!, 0.138 mins elapsed.
## 2020-12-03 18:26:05 : Adding s1Inj to PeakMatrix for Chr (5 of 20)!, 0.18 mins elapsed.
## 2020-12-03 18:26:07 : Adding s1Inj to PeakMatrix for Chr (6 of 20)!, 0.222 mins elapsed.
## 2020-12-03 18:26:10 : Adding s1Inj to PeakMatrix for Chr (7 of 20)!, 0.261 mins elapsed.
## 2020-12-03 18:26:12 : Adding s1Inj to PeakMatrix for Chr (8 of 20)!, 0.304 mins elapsed.
## 2020-12-03 18:26:15 : Adding s1Inj to PeakMatrix for Chr (9 of 20)!, 0.344 mins elapsed.
## 2020-12-03 18:26:17 : Adding s1Inj to PeakMatrix for Chr (10 of 20)!, 0.384 mins elapsed.
## 2020-12-03 18:26:20 : Adding s1Inj to PeakMatrix for Chr (11 of 20)!, 0.423 mins elapsed.
## 2020-12-03 18:26:22 : Adding s1Inj to PeakMatrix for Chr (12 of 20)!, 0.466 mins elapsed.
## 2020-12-03 18:26:24 : Adding s1Inj to PeakMatrix for Chr (13 of 20)!, 0.503 mins elapsed.
## 2020-12-03 18:26:27 : Adding s1Inj to PeakMatrix for Chr (14 of 20)!, 0.541 mins elapsed.
## 2020-12-03 18:26:29 : Adding s1Inj to PeakMatrix for Chr (15 of 20)!, 0.578 mins elapsed.
## 2020-12-03 18:26:31 : Adding s1Inj to PeakMatrix for Chr (16 of 20)!, 0.616 mins elapsed.
## 2020-12-03 18:26:33 : Adding s1Inj to PeakMatrix for Chr (17 of 20)!, 0.652 mins elapsed.
## 2020-12-03 18:26:36 : Adding s1Inj to PeakMatrix for Chr (18 of 20)!, 0.69 mins elapsed.
## 2020-12-03 18:26:38 : Adding s1Inj to PeakMatrix for Chr (19 of 20)!, 0.726 mins elapsed.
## 2020-12-03 18:26:40 : Adding s1Inj to PeakMatrix for Chr (20 of 20)!, 0.761 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2020-12-03 18:26:49 : Adding s3Inj to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2020-12-03 18:26:51 : Adding s3Inj to PeakMatrix for Chr (2 of 20)!, 0.047 mins elapsed.
## 2020-12-03 18:26:54 : Adding s3Inj to PeakMatrix for Chr (3 of 20)!, 0.085 mins elapsed.
## 2020-12-03 18:26:56 : Adding s3Inj to PeakMatrix for Chr (4 of 20)!, 0.12 mins elapsed.
## 2020-12-03 18:26:58 : Adding s3Inj to PeakMatrix for Chr (5 of 20)!, 0.158 mins elapsed.
## 2020-12-03 18:27:00 : Adding s3Inj to PeakMatrix for Chr (6 of 20)!, 0.195 mins elapsed.
## 2020-12-03 18:27:02 : Adding s3Inj to PeakMatrix for Chr (7 of 20)!, 0.23 mins elapsed.
## 2020-12-03 18:27:05 : Adding s3Inj to PeakMatrix for Chr (8 of 20)!, 0.268 mins elapsed.
## 2020-12-03 18:27:07 : Adding s3Inj to PeakMatrix for Chr (9 of 20)!, 0.304 mins elapsed.
## 2020-12-03 18:27:09 : Adding s3Inj to PeakMatrix for Chr (10 of 20)!, 0.339 mins elapsed.
## 2020-12-03 18:27:11 : Adding s3Inj to PeakMatrix for Chr (11 of 20)!, 0.375 mins elapsed.
## 2020-12-03 18:27:13 : Adding s3Inj to PeakMatrix for Chr (12 of 20)!, 0.412 mins elapsed.
## 2020-12-03 18:27:15 : Adding s3Inj to PeakMatrix for Chr (13 of 20)!, 0.446 mins elapsed.
## 2020-12-03 18:27:17 : Adding s3Inj to PeakMatrix for Chr (14 of 20)!, 0.481 mins elapsed.
## 2020-12-03 18:27:20 : Adding s3Inj to PeakMatrix for Chr (15 of 20)!, 0.515 mins elapsed.
## 2020-12-03 18:27:22 : Adding s3Inj to PeakMatrix for Chr (16 of 20)!, 0.55 mins elapsed.
## 2020-12-03 18:27:24 : Adding s3Inj to PeakMatrix for Chr (17 of 20)!, 0.584 mins elapsed.
## 2020-12-03 18:27:26 : Adding s3Inj to PeakMatrix for Chr (18 of 20)!, 0.618 mins elapsed.
## 2020-12-03 18:27:28 : Adding s3Inj to PeakMatrix for Chr (19 of 20)!, 0.651 mins elapsed.
## 2020-12-03 18:27:30 : Adding s3Inj to PeakMatrix for Chr (20 of 20)!, 0.684 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2020-12-03 18:27:38 : Adding s3UI to PeakMatrix for Chr (1 of 20)!, 0.009 mins elapsed.
## 2020-12-03 18:27:41 : Adding s3UI to PeakMatrix for Chr (2 of 20)!, 0.047 mins elapsed.
## 2020-12-03 18:27:43 : Adding s3UI to PeakMatrix for Chr (3 of 20)!, 0.085 mins elapsed.
## 2020-12-03 18:27:45 : Adding s3UI to PeakMatrix for Chr (4 of 20)!, 0.12 mins elapsed.
## 2020-12-03 18:27:47 : Adding s3UI to PeakMatrix for Chr (5 of 20)!, 0.158 mins elapsed.
## 2020-12-03 18:27:49 : Adding s3UI to PeakMatrix for Chr (6 of 20)!, 0.195 mins elapsed.
## 2020-12-03 18:27:52 : Adding s3UI to PeakMatrix for Chr (7 of 20)!, 0.231 mins elapsed.
## 2020-12-03 18:27:54 : Adding s3UI to PeakMatrix for Chr (8 of 20)!, 0.269 mins elapsed.
## 2020-12-03 18:27:56 : Adding s3UI to PeakMatrix for Chr (9 of 20)!, 0.306 mins elapsed.
## 2020-12-03 18:27:58 : Adding s3UI to PeakMatrix for Chr (10 of 20)!, 0.342 mins elapsed.
## 2020-12-03 18:28:00 : Adding s3UI to PeakMatrix for Chr (11 of 20)!, 0.378 mins elapsed.
## 2020-12-03 18:28:03 : Adding s3UI to PeakMatrix for Chr (12 of 20)!, 0.416 mins elapsed.
## 2020-12-03 18:28:05 : Adding s3UI to PeakMatrix for Chr (13 of 20)!, 0.451 mins elapsed.
## 2020-12-03 18:28:07 : Adding s3UI to PeakMatrix for Chr (14 of 20)!, 0.485 mins elapsed.
## 2020-12-03 18:28:09 : Adding s3UI to PeakMatrix for Chr (15 of 20)!, 0.519 mins elapsed.
## 2020-12-03 18:28:11 : Adding s3UI to PeakMatrix for Chr (16 of 20)!, 0.554 mins elapsed.
## 2020-12-03 18:28:13 : Adding s3UI to PeakMatrix for Chr (17 of 20)!, 0.587 mins elapsed.
## 2020-12-03 18:28:15 : Adding s3UI to PeakMatrix for Chr (18 of 20)!, 0.622 mins elapsed.
## 2020-12-03 18:28:17 : Adding s3UI to PeakMatrix for Chr (19 of 20)!, 0.656 mins elapsed.
## 2020-12-03 18:28:19 : Adding s3UI to PeakMatrix for Chr (20 of 20)!, 0.689 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2020-12-03 18:28:28 : Adding s2UI to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2020-12-03 18:28:30 : Adding s2UI to PeakMatrix for Chr (2 of 20)!, 0.049 mins elapsed.
## 2020-12-03 18:28:33 : Adding s2UI to PeakMatrix for Chr (3 of 20)!, 0.089 mins elapsed.
## 2020-12-03 18:28:35 : Adding s2UI to PeakMatrix for Chr (4 of 20)!, 0.127 mins elapsed.
## 2020-12-03 18:28:37 : Adding s2UI to PeakMatrix for Chr (5 of 20)!, 0.165 mins elapsed.
## 2020-12-03 18:28:40 : Adding s2UI to PeakMatrix for Chr (6 of 20)!, 0.204 mins elapsed.
## 2020-12-03 18:28:42 : Adding s2UI to PeakMatrix for Chr (7 of 20)!, 0.241 mins elapsed.
## 2020-12-03 18:28:44 : Adding s2UI to PeakMatrix for Chr (8 of 20)!, 0.283 mins elapsed.
## 2020-12-03 18:28:47 : Adding s2UI to PeakMatrix for Chr (9 of 20)!, 0.334 mins elapsed.
## 2020-12-03 18:28:50 : Adding s2UI to PeakMatrix for Chr (10 of 20)!, 0.382 mins elapsed.
## 2020-12-03 18:28:53 : Adding s2UI to PeakMatrix for Chr (11 of 20)!, 0.431 mins elapsed.
## 2020-12-03 18:28:56 : Adding s2UI to PeakMatrix for Chr (12 of 20)!, 0.47 mins elapsed.
## 2020-12-03 18:28:58 : Adding s2UI to PeakMatrix for Chr (13 of 20)!, 0.505 mins elapsed.
## 2020-12-03 18:29:00 : Adding s2UI to PeakMatrix for Chr (14 of 20)!, 0.544 mins elapsed.
## 2020-12-03 18:29:02 : Adding s2UI to PeakMatrix for Chr (15 of 20)!, 0.579 mins elapsed.
## 2020-12-03 18:29:04 : Adding s2UI to PeakMatrix for Chr (16 of 20)!, 0.613 mins elapsed.
## 2020-12-03 18:29:07 : Adding s2UI to PeakMatrix for Chr (17 of 20)!, 0.653 mins elapsed.
## 2020-12-03 18:29:09 : Adding s2UI to PeakMatrix for Chr (18 of 20)!, 0.689 mins elapsed.
## 2020-12-03 18:29:11 : Adding s2UI to PeakMatrix for Chr (19 of 20)!, 0.731 mins elapsed.
## 2020-12-03 18:29:14 : Adding s2UI to PeakMatrix for Chr (20 of 20)!, 0.772 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2020-12-03 18:29:27 : Adding s2Inj to PeakMatrix for Chr (1 of 20)!, 0.011 mins elapsed.
## 2020-12-03 18:29:29 : Adding s2Inj to PeakMatrix for Chr (2 of 20)!, 0.048 mins elapsed.
## 2020-12-03 18:29:31 : Adding s2Inj to PeakMatrix for Chr (3 of 20)!, 0.086 mins elapsed.
## 2020-12-03 18:29:33 : Adding s2Inj to PeakMatrix for Chr (4 of 20)!, 0.124 mins elapsed.
## 2020-12-03 18:29:36 : Adding s2Inj to PeakMatrix for Chr (5 of 20)!, 0.16 mins elapsed.
## 2020-12-03 18:29:38 : Adding s2Inj to PeakMatrix for Chr (6 of 20)!, 0.195 mins elapsed.
## 2020-12-03 18:29:40 : Adding s2Inj to PeakMatrix for Chr (7 of 20)!, 0.23 mins elapsed.
## 2020-12-03 18:29:42 : Adding s2Inj to PeakMatrix for Chr (8 of 20)!, 0.268 mins elapsed.
## 2020-12-03 18:29:44 : Adding s2Inj to PeakMatrix for Chr (9 of 20)!, 0.303 mins elapsed.
## 2020-12-03 18:29:46 : Adding s2Inj to PeakMatrix for Chr (10 of 20)!, 0.338 mins elapsed.
## 2020-12-03 18:29:48 : Adding s2Inj to PeakMatrix for Chr (11 of 20)!, 0.373 mins elapsed.
## 2020-12-03 18:29:51 : Adding s2Inj to PeakMatrix for Chr (12 of 20)!, 0.409 mins elapsed.
## 2020-12-03 18:29:53 : Adding s2Inj to PeakMatrix for Chr (13 of 20)!, 0.442 mins elapsed.
## 2020-12-03 18:29:55 : Adding s2Inj to PeakMatrix for Chr (14 of 20)!, 0.476 mins elapsed.
## 2020-12-03 18:29:57 : Adding s2Inj to PeakMatrix for Chr (15 of 20)!, 0.509 mins elapsed.
## 2020-12-03 18:29:59 : Adding s2Inj to PeakMatrix for Chr (16 of 20)!, 0.543 mins elapsed.
## 2020-12-03 18:30:01 : Adding s2Inj to PeakMatrix for Chr (17 of 20)!, 0.576 mins elapsed.
## 2020-12-03 18:30:03 : Adding s2Inj to PeakMatrix for Chr (18 of 20)!, 0.61 mins elapsed.
## 2020-12-03 18:30:05 : Adding s2Inj to PeakMatrix for Chr (19 of 20)!, 0.643 mins elapsed.
## 2020-12-03 18:30:07 : Adding s2Inj to PeakMatrix for Chr (20 of 20)!, 0.675 mins elapsed.
## .createArrowGroup : Arrow Group already exists! Dropping Group from ArrowFile! This will take ~10-30 seconds!
## .dropGroupsFromArrow : Initializing Temp ArrowFile
## .dropGroupsFromArrow : Adding Metadata to Temp ArrowFile
## .dropGroupsFromArrow : Adding SubGroups to Temp ArrowFile
## .dropGroupsFromArrow : Move Temp ArrowFile to ArrowFile
## 2020-12-03 18:30:15 : Adding s1UI to PeakMatrix for Chr (1 of 20)!, 0.01 mins elapsed.
## 2020-12-03 18:30:17 : Adding s1UI to PeakMatrix for Chr (2 of 20)!, 0.044 mins elapsed.
## 2020-12-03 18:30:19 : Adding s1UI to PeakMatrix for Chr (3 of 20)!, 0.078 mins elapsed.
## 2020-12-03 18:30:21 : Adding s1UI to PeakMatrix for Chr (4 of 20)!, 0.111 mins elapsed.
## 2020-12-03 18:30:24 : Adding s1UI to PeakMatrix for Chr (5 of 20)!, 0.145 mins elapsed.
## 2020-12-03 18:30:26 : Adding s1UI to PeakMatrix for Chr (6 of 20)!, 0.179 mins elapsed.
## 2020-12-03 18:30:28 : Adding s1UI to PeakMatrix for Chr (7 of 20)!, 0.212 mins elapsed.
## 2020-12-03 18:30:30 : Adding s1UI to PeakMatrix for Chr (8 of 20)!, 0.246 mins elapsed.
## 2020-12-03 18:30:32 : Adding s1UI to PeakMatrix for Chr (9 of 20)!, 0.279 mins elapsed.
## 2020-12-03 18:30:34 : Adding s1UI to PeakMatrix for Chr (10 of 20)!, 0.312 mins elapsed.
## 2020-12-03 18:30:36 : Adding s1UI to PeakMatrix for Chr (11 of 20)!, 0.346 mins elapsed.
## 2020-12-03 18:30:38 : Adding s1UI to PeakMatrix for Chr (12 of 20)!, 0.384 mins elapsed.
## 2020-12-03 18:30:40 : Adding s1UI to PeakMatrix for Chr (13 of 20)!, 0.417 mins elapsed.
## 2020-12-03 18:30:42 : Adding s1UI to PeakMatrix for Chr (14 of 20)!, 0.45 mins elapsed.
## 2020-12-03 18:30:44 : Adding s1UI to PeakMatrix for Chr (15 of 20)!, 0.482 mins elapsed.
## 2020-12-03 18:30:46 : Adding s1UI to PeakMatrix for Chr (16 of 20)!, 0.515 mins elapsed.
## 2020-12-03 18:30:48 : Adding s1UI to PeakMatrix for Chr (17 of 20)!, 0.549 mins elapsed.
## 2020-12-03 18:30:50 : Adding s1UI to PeakMatrix for Chr (18 of 20)!, 0.583 mins elapsed.
## 2020-12-03 18:30:52 : Adding s1UI to PeakMatrix for Chr (19 of 20)!, 0.616 mins elapsed.
## 2020-12-03 18:30:54 : Adding s1UI to PeakMatrix for Chr (20 of 20)!, 0.647 mins elapsed.
## Overriding previous entry for ReadsInPeaks
## Overriding previous entry for FRIP
## ArchR logging successful to : ArchRLogs/ArchR-addPeakMatrix-175b85a2b745d-Date-2020-12-03_Time-18-25-47.log## get ArchR markers for motif enrichment
markerPeaks_subHybrid <- getMarkerFeatures(
ArchRProj = archProj,
useMatrix = "PeakMatrix",
groupBy = "predAnn",
bias = c("TSSEnrichment", "log10(nFrags)"),
testMethod = "wilcoxon")## ArchR logging to : ArchRLogs/ArchR-getMarkerFeatures-175b82744b83f-Date-2020-12-03_Time-18-30-56.log
## If there is an issue, please report to github with logFile!
## MatrixClass = Sparse.Integer.Matrix
## 2020-12-03 18:30:56 : Matching Known Biases, 0.004 mins elapsed.
## 2020-12-03 18:30:58 : Computing Pairwise Tests (1 of 5), 0.031 mins elapsed.
## 2020-12-03 18:31:26 : Computing Pairwise Tests (2 of 5), 0.494 mins elapsed.
## 2020-12-03 18:31:53 : Computing Pairwise Tests (3 of 5), 0.95 mins elapsed.
## 2020-12-03 18:32:18 : Computing Pairwise Tests (4 of 5), 1.374 mins elapsed.
## 2020-12-03 18:32:44 : Computing Pairwise Tests (5 of 5), 1.804 mins elapsed.
## ###########
## 2020-12-03 18:33:11 : Completed Pairwise Tests, 2.255 mins elapsed.
## ###########
## ArchR logging successful to : ArchRLogs/ArchR-getMarkerFeatures-175b82744b83f-Date-2020-12-03_Time-18-30-56.logpeakMarkers_subHybrid <- getMarkers(markerPeaks_subHybrid,
cutOff = "FDR <= 0.01 & Log2FC >= 1.25",
returnGR = TRUE)
heatmapMarkerPeaks <- plotMarkerHeatmap(
seMarker = markerPeaks_subHybrid,
cutOff = "FDR <= 0.01 & Log2FC >= 1.25",
transpose = TRUE
)## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-175b8109c443d-Date-2020-12-03_Time-18-33-13.log
## If there is an issue, please report to github with logFile!
## Identified 20658 markers!
## Preparing Main Heatmap..
## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-175b8109c443d-Date-2020-12-03_Time-18-33-13.logComplexHeatmap::draw(heatmapMarkerPeaks, heatmap_legend_side = "bot", annotation_legend_side = "bot")
Color RNA-seq data based on new markers
plotByExpr2 <- function(rna, genes, counts, name="feature", scale=TRUE){
genes <- genes[genes %in% rownames(countHBC)]
if(scale){
scaledCounts <- t(scale(t(counts[genes,]+1e-5)))
scaledCounts <- scaledCounts[!is.na(rowVars(scaledCounts)),]
yhat <- colSums(scaledCounts)
} else {
yhat <- colSums(counts[genes,]+1e-5)
}
rna <- AddMetaData(rna, yhat, col.name = name)
# featureplot
p1 <- FeaturePlot(rna, features = name) + theme(legend.position = "none")
# boxplot
df <- data.frame(y=yhat, cluster=rna$ctHybridSub)
p2 <- ggplot(df, aes(x=cluster, y=yhat)) +
geom_boxplot() +
theme_classic()
cowplot::plot_grid(p1, p2, nrow=1)
}
peakSet <- archProj@peakSet
markersHybrid <- peakMarkers_subHybrid$C2_Hybrid
markersHybrid1 <- peakMarkers_subHybrid$hybrid1
markersHybrid2 <- peakMarkers_subHybrid$hybrid2
peaksHybrid <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid, type = "equal"))]
stopifnot(length(peaksHybrid) == length(markersHybrid))
peaksHybridFunnel <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid1, type = "equal"))]
stopifnot(length(peaksHybridFunnel) == length(markersHybrid1))
peaksHybridHorn <- peakSet[queryHits(findOverlaps(peakSet, markersHybrid2, type = "equal"))]
stopifnot(length(peaksHybridHorn) == length(markersHybrid2))
plotByExpr2(rna,
genes=intersect(unique(mcols(peaksHybrid)$nearestGene), rownames(countHBC)),
counts=countHBC,
name="HybridMarkers",
scale=TRUE)
plotByExpr2(rna,
genes=intersect(unique(mcols(peaksHybridFunnel)$nearestGene), rownames(countHBC)),
counts=countHBC,
name="Hybrid1Markers",
scale=TRUE)
plotByExpr2(rna,
genes=intersect(unique(mcols(peaksHybridHorn)$nearestGene), rownames(countHBC)),
counts=countHBC,
name="Hybrid2Markers",
scale=TRUE)
Motif enrichment for sub-hyrbid clusters
archProj <- addMotifAnnotations(ArchRProj = archProj,
motifSet = "cisbp",
name = "Motif",
force=TRUE)## No methods found in package 'IRanges' for request: 'score' when loading 'TFBSTools'## ArchR logging to : ArchRLogs/ArchR-addMotifAnnotations-175b825955b7a-Date-2020-12-03_Time-18-33-25.log
## If there is an issue, please report to github with logFile!## peakAnnotation name already exists! Overriding.## 2020-12-03 18:33:26 : Gettting Motif Set, Species : Mus musculus, 0.002 mins elapsed.## Using version 2 motifs!## 2020-12-03 18:33:29 : Finding Motif Positions with motifmatchr!, 0.039 mins elapsed.## 2020-12-03 18:36:37 : Creating Motif Overlap Matrix, 3.175 mins elapsed.## 2020-12-03 18:36:40 : Finished Getting Motif Info!, 3.224 mins elapsed.## ArchR logging successful to : ArchRLogs/ArchR-addMotifAnnotations-175b825955b7a-Date-2020-12-03_Time-18-33-25.log## upregulated motifs
motifsUp <- peakAnnoEnrichment(
seMarker = markerPeaks_subHybrid,
ArchRProj = archProj,
peakAnnotation = "Motif",
cutOff = "FDR <= 0.05 & Log2FC >= 0.5"
)## ArchR logging to : ArchRLogs/ArchR-peakAnnoEnrichment-175b870a8b9dc-Date-2020-12-03_Time-18-36-43.log
## If there is an issue, please report to github with logFile!## 2020-12-03 18:36:48 : Computing Enrichments 1 of 5, 0.075 mins elapsed.## 2020-12-03 18:36:48 : Computing Enrichments 2 of 5, 0.084 mins elapsed.## 2020-12-03 18:36:49 : Computing Enrichments 3 of 5, 0.094 mins elapsed.## 2020-12-03 18:36:50 : Computing Enrichments 4 of 5, 0.103 mins elapsed.## 2020-12-03 18:36:50 : Computing Enrichments 5 of 5, 0.113 mins elapsed.## ArchR logging successful to : ArchRLogs/ArchR-peakAnnoEnrichment-175b870a8b9dc-Date-2020-12-03_Time-18-36-43.log# the top 10 TFs are same for all 3 hybrid clusters which is why there's only 10 cols for all 3 clusters.
heatmapMotifsUp <- plotEnrichHeatmap(motifsUp, transpose = TRUE)## ArchR logging to : ArchRLogs/ArchR-plotEnrichHeatmap-175b85542544c-Date-2020-12-03_Time-18-36-51.log
## If there is an issue, please report to github with logFile!## Adding Annotations..## Preparing Main Heatmap..ComplexHeatmap::draw(heatmapMotifsUp, heatmap_legend_side = "bot", annotation_legend_side = "bot")
# note that Gm5294 is an ortholog of the FOXL3-OT1 lncRNA and paralog of FOXL1
getTopDf <- function(motifsRes, cl){
df <- data.frame(TF = rownames(motifsRes), mlog10Padj = assay(motifsRes)[,cl])
df <- df[order(df$mlog10Padj, decreasing = TRUE),]
df$rank <- seq_len(nrow(df))
return(df)
}
plotMotifs <- function(motifsRes, cl){
df <- getTopDf(motifsRes, cl)
ggUp <- ggplot(df, aes(rank, mlog10Padj, color = mlog10Padj)) +
geom_point(size = 1) +
ggrepel::geom_label_repel(
data = df[rev(seq_len(30)), ], aes(x = rank, y = mlog10Padj, label = TF),
size = 1.5,
nudge_x = 2,
color = "black"
) + theme_ArchR() +
ylab("-log10(P-adj) Motif Enrichment") +
xlab("Rank Sorted TFs Enriched") +
scale_color_gradientn(colors = paletteContinuous(set = "comet"))
ggUp
}
# the upregulated markers for all hybrid clusters are enriched in Fox motifs.
plotMotifs(motifsUp, "C2_Hybrid") + ggtitle("Hybrid")
plotMotifs(motifsUp, "hybrid1") + ggtitle("hybrid1")
plotMotifs(motifsUp, "hybrid1") + ggtitle("hybrid1")
# Fox motifs are not enriched in either resting or activated cell types.
plotMotifs(motifsUp, "resting") + ggtitle("Resting")
plotMotifs(motifsUp, "C1_Inj") + ggtitle("Activated")
Smarcc1 motifs
Smarcc1 was recently found to regulate Keratin genes and thereby cellular fate in mammalian embryogenesis: https://www.nature.com/articles/s41586-020-2647-4. Here, we find that the markers for the activated cells are most enriched in Smarcc1 motifs. To dig a little deeper, we therefore check which of these marker peaks that may be close to Krt genes have Smarcc1 motifs here.
#detach("package:matrixStats", unload=TRUE)
# get the marker peaks for activated cells
actMarkers <- peakMarkers_subHybrid$C1_Inj
# subset to those that have Smarcc1 motif
mot <- archProj@peakAnnotation$Motif
matches <- readRDS(mot$Matches)
markerHits <- queryHits(findOverlaps(query = SummarizedExperiment::rowRanges(matches),
subject = actMarkers,
type = "equal"))
smarcHits <- which(assays(matches, "matches")[[1]][,grep(x=colnames(assays(matches, "matches")[[1]]), pattern="Smarcc1")])
markerSmarcRanges <- SummarizedExperiment::rowRanges(matches)[intersect(markerHits, smarcHits)]
# check which Krt genes are in there
krtWithMotif <- sort(unname(grep(x=mcols(markerSmarcRanges)$nearestGene, pattern="Krt", value=TRUE)))
krtWithMotif## [1] "Krt13" "Krt14" "Krt16" "Krt18" "Krt19" "Krt20"
## [7] "Krt20" "Krt26" "Krt5" "Krt5" "Krt6a" "Krt84"
## [13] "Krtap11-1" "Krtap2-4" "Krtap4-16"# plot their accessibility and expression
for(kk in 1:length(krtWithMotif)){
gene <- krtWithMotif[kk]
if(gene %in% rownames(rna)){
prna <- FeaturePlot(rna, features=gene)
patac <- FeaturePlot(atac, features=gene) + xlim(c(-6,6))
print(cowplot::plot_grid(prna, patac, nrow=1))
} else {
print(FeaturePlot(atac, features=gene) + xlim(c(-6,6)))
}
}## Warning: Could not find Krt13 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krt14 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt16 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt18 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt19 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt20 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt20 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt26 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt6a in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 55 rows containing missing values (geom_point).
## Warning: Could not find Krt84 in the default search locations, found in ACTIVITY
## assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krtap11-1 in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krtap2-4 in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).## Warning: Could not find Krtap4-16 in the default search locations, found in
## ACTIVITY assay instead## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 55 rows containing missing values (geom_point).
Check expression of Fox TFs and target genes in scRNA-seq data
Peaks that are more accessible in hybrid cells are upregulated in motifs bound by Fox transcription factors. We expect Fox transcription factors to be higher expressed in hybrid clusters. Since Fox transcription factors generally are repressors, we could expect the target genes of Fox transcription factors to be lower expressed in the hybrid clusters.
Evidence from the literature: - Foxg1 expression is in horizontal basal cells of OE, Foxg1 -/- mice are born with defects in the OE: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2878634/ - Fox family TF motifs are most enriched in olfactory receptors: https://arxiv.org/pdf/1201.2933.pdf - Foxl1 are markers of microvillous cells: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2878634/
# selected Fox TFs as identified with motif enrichment
# top 10 motifs TFs from hybrid cells:
selectedFoxTfs <- c("Floxl1", "Foxa1", "Foxb2", "Foxb1", "Foxc1", "Foxc2", "Foxs1", "Foxd1")
selectedFoxTfs <- intersect(selectedFoxTfs, rownames(countHBC))
countsFoxSelected <- countHBC[selectedFoxTfs,]
foxSumSelected <- colSums(countsFoxSelected)
rna <- AddMetaData(rna, foxSumSelected, col.name = "foxSumSelected")
pFoxSum <- FeaturePlot(rna, features = "foxSumSelected")
df <- data.frame(cl=rna$ctHybridSub,
foxSum=foxSumSelected)
pBoxFoxSum <- ggplot(df, aes(x=cl, y=foxSum)) +
geom_boxplot() +
theme_classic() +
ylab("Sum of selected Fox TFs expression")
cowplot::plot_grid(pFoxSum, pBoxFoxSum)
#ggsave("../plots/foxSumRNA.png", width=8)
## scaled
scaledCountsFoxSelected <- t(scale(t(countHBC[selectedFoxTfs,]+1e-5)))
scaledCountsFoxSelected <- scaledCountsFoxSelected[!is.na(rowVars(scaledCountsFoxSelected)),]
foxSumScaledSelected <- colSums(scaledCountsFoxSelected)
rna <- AddMetaData(rna, foxSumScaledSelected, col.name = "foxSumSelected_scaled")
pFoxSumScaled <- FeaturePlot(rna, features = "foxSumSelected_scaled")
df <- data.frame(cl=rna$ctHybridSub,
foxScaled=foxSumScaledSelected)
pBoxFoxSumScaled <- ggplot(df, aes(x=cl, y=foxScaled)) +
geom_boxplot() +
theme_classic() +
ylab("Sum of selected Fox TFs scaled expression")
cowplot::plot_grid(pFoxSumScaled, pBoxFoxSumScaled)
# ggsave("../plots/foxSumRNAScaled.png", width=8)Composite plot
## DR of RNA
pRNA1_fin <- pRNA1 + ggtitle("scRNA-seq data") +
theme(legend.title = element_text(size = 12),
legend.text = element_text(size = 10),
legend.key.size = unit(1/2, "cm"))
## DR of RNA with ATAC markers expression highlighted
pMarkersATAC <- cowplot::plot_grid(plotlist=list(pc1,pc2,pc3), nrow=1)
pMarkersATAC
## RNA and ATAC cell type transfer
pCTrans <- cowplot::plot_grid(ppred4_2 + NoLegend(), ppred2_2, nrow=1)
## markers and motifs for subhybrid clusters
## marker
heatmapMarkerPeaksMatrix <- plotMarkerHeatmap(
seMarker = markerPeaks_subHybrid,
cutOff = "FDR <= 0.01 & Log2FC >= 1.25",
transpose = TRUE,
returnMat = TRUE
)## ArchR logging to : ArchRLogs/ArchR-plotMarkerHeatmap-175b835e97d00-Date-2020-12-03_Time-18-37-30.log
## If there is an issue, please report to github with logFile!## Identified 20658 markers!## ArchR logging successful to : ArchRLogs/ArchR-plotMarkerHeatmap-175b835e97d00-Date-2020-12-03_Time-18-37-30.logpMarker <- Heatmap(heatmapMarkerPeaksMatrix,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = paletteContinuous(set = "solarExtra", n = 100),
show_column_names = FALSE,
row_names_gp = gpar(fontsize = 13),
heatmap_legend_param = list(title = "Z-score"))
## motif
heatmapMotifsMatrix <- plotEnrichHeatmap(motifsUp, transpose = TRUE, returnMatrix = TRUE)## ArchR logging to : ArchRLogs/ArchR-plotEnrichHeatmap-175b8740da6a6-Date-2020-12-03_Time-18-37-32.log
## If there is an issue, please report to github with logFile!colnames(heatmapMotifsMatrix) <- unlist(lapply(strsplit(colnames(heatmapMotifsMatrix), split="_"), "[[", 1))
pMotif <- Heatmap(heatmapMotifsMatrix,
cluster_rows = FALSE,
cluster_columns = FALSE,
col = paletteContinuous(set = "comet", n = 100),
column_names_gp = gpar(fontsize = 11),
row_names_gp = gpar(fontsize = 13),
heatmap_legend_param = list(title = "-log10\n p-adj."))
pHeat <- cowplot::plot_grid(grid::grid.grabExpr(ComplexHeatmap::draw(pMarker)),
grid::grid.grabExpr(ComplexHeatmap::draw(pMotif)),
nrow = 2)
pComp1 <- cowplot::plot_grid(pRNA1_fin, pMarkersATAC, nrow=1,
rel_widths = c(0.4, 1), labels = c("a", "b"))
pComp2 <- cowplot::plot_grid(pCTrans, pHeat, ncol = 1,
rel_heights = c(0.8, 1), labels = c("c", "d"))
cowplot::plot_grid(pComp1, pComp2, nrow = 2,
rel_widths= c(1/3, 1))
# ggsave("../plots/integrationComposite.pdf", height=13, width=12)
# ggsave("../plots/integrationComposite.png", height=13, width=12)old plots
Explore the activated HBC RNA-seq clusters
Macrophage markers
FeaturePlot(atac, features = "Msr1")## Warning: Could not find Msr1 in the default search locations, found in ACTIVITY
## assay instead
FeaturePlot(rna, features = "Msr1")
Cell cycle genes
download.file("https://raw.githubusercontent.com/rufletch/p63-HBC-diff/master/ref/cellCycle40.txt", destfile = "~/tmp/ccGenes40.txt")
cc40Genes <- read.table("~/tmp/ccGenes40.txt", stringsAsFactors = FALSE)[,1]
countsCC40 <- countHBC[cc40Genes,]+1e-5
cc40Sum <- colSums(countsCC40)
rna <- AddMetaData(rna, cc40Sum, col.name = "cc40Sum")
FeaturePlot(rna, features = "cc40Sum")
## scaled
scaledCountsCC40 <- t(scale(t(countHBC[cc40Genes,]+1e-5)))
scaledCountsCC40 <- scaledCountsCC40[!is.na(rowVars(scaledCountsCC40)),]
cc40Sum <- colSums(scaledCountsCC40)
rna <- AddMetaData(rna, cc40Sum, col.name = "cc40Sum_scaled")
FeaturePlot(rna, features = "cc40Sum_scaled")
download.file("https://raw.githubusercontent.com/rufletch/p63-HBC-diff/master/ref/cell_cycle.txt", destfile = "~/tmp/ccGenes.txt")
ccGenes <- read.table("~/tmp/ccGenes.txt", stringsAsFactors = FALSE)[,1]
ccGenes <- ccGenes[ccGenes %in% rownames(countHBC)]
countsCC <- countHBC[ccGenes,]+1e-5
ccSum <- colSums(countsCC)
rna <- AddMetaData(rna, ccSum, col.name = "ccSum")
FeaturePlot(rna, features = "ccSum")
## scaled
scaledCountsCC <- t(scale(t(countHBC[ccGenes,]+1e-5)))
scaledCountsCC <- scaledCountsCC[!is.na(rowVars(scaledCountsCC)),]
ccSum <- colSums(scaledCountsCC)
rna <- AddMetaData(rna, ccSum, col.name = "ccSum_scaled")
FeaturePlot(rna, features = "ccSum_scaled")
Wound response genes
Genes obtained from paper on HBC injury, Fig 3D caption: `Krt6a, Krt16, Sprr1a, Sprr2a3, Krtdap, Dmkn, Sbsn, and Hbegf’
wrGenes <- c("Krt6a", "Krt16", "Sprr1a", "Sprr2a3", "Krtdap", "Dmkn", "Sbsn", "Hbegf")
countsWR <- countHBC[wrGenes,]+1e-5
wrSum <- colSums(countsWR)
rna <- AddMetaData(rna, wrSum, col.name = "wrSum")
FeaturePlot(rna, features = "wrSum")
## scaled
scaledCountsWR <- t(scale(t(countHBC[wrGenes,]+1e-5)))
scaledCountsWR <- scaledCountsWR[!is.na(rowVars(scaledCountsWR)),]
wrSum <- colSums(scaledCountsWR)
rna <- AddMetaData(rna, wrSum, col.name = "wrSum_scaled")
FeaturePlot(rna, features = "wrSum_scaled")
Respiratory markers
respGenesRebecca <- c("Gp2", "Cyp4a12b", "Kcnj16", "Tff2", "Chad", "Acsm1", "Kcne3", "Vtcn1", "Kcnj15", "Cldn8", "Pigr", "Muc5b", "Calml3", "Kng2", "Cyp2b10")
## note that Cxcl14 and Meg3 are only described for human samples in the Methods.
respGenesDavid <- c("Foxj1", "Muc5ac", "Cxcl14", "Meg3")
## genes used by Boying
respGenesBoying <- c("Reg3g", "Muc5ac", "Muc5b", "Meg3", "Krt5", "Krt14", "Trp63")
FeaturePlot(rna, features = respGenesRebecca[1:9])
FeaturePlot(rna, features = respGenesRebecca[10:15])
FeaturePlot(rna, features = respGenesDavid)
FeaturePlot(rna, features = respGenesBoying) + NoLegend()